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Mrk 421 as a case study for TeV and X-ray variability in 
leptohadronic models 



m 

o 

o 



Oh! 

6 



> 

in 

o 
m 



1 



A. Mastichiadis (AM) * M. Petropoulou (MP) f and S. Dimitrakoudis (SD) J 

Department of Physics, University of Athens, Panepistimiopolis, GR 15783 Zografos, Greece 



Received. . ./Accepted. . 



ABSTRACT 

We investigate the origin of high-energy emission in blazars within the context of 
the leptohadronic one-zone model. We find that 7-ray emission can be attributed to 
synchrotron radiation either from protons or from secondary leptons produced via 
photohadronic processes. These possibilities imply differences not only in the spectral 
energy distribution (SED) but also in the variability signatures, especially in the X- 
and 7-ray regime. Thus, the temporal behavior of each leptohadronic scenario can 
be used to probe the particle population responsible for the high-energy emission as 
it can give extra information not available by spectral fits. In the present work we 
apply these ideas to the non-thermal emission of Mrk 421, which is one of the best 
monitored TeV blazars. We focus on the observations of March 2001, since during that 
period Mrk 421 showed multiple flares that have been observed in detail both in X- 
rays and 7-rays. First, we obtain pre-flaring fits to the SED using the different types 
of leptohadronic scenarios. Then, we introduce random-walk type, small-amplitude 
variations on the injection compactness or on the maximum energy of radiating par- 
ticles and follow the subsequent response of the radiated photon spectrum. For each 
leptohadronic scenario, we calculate the X-ray and 7-ray fluxes and investigate their 
possible correlation. Whenever the 'input' variations lead, apart from flux variability, 
also to spectral variability, we present the resulting relations between the spectral in- 
dex and the flux, both in X-rays and 7-rays. We find that proton synchrotron models 
are favoured energetically but require fine tuning between electron and proton param- 
eters to reproduce the observed quadratic behaviour between X-rays and TeV 7— rays. 
On the other hand, models based on pion-decay can reproduce this behaviour in a 
much more natural way. 

Key words: astroparticle physics - radiation mechanisms: non-thermal - gamma 
rays: galaxies - galaxies: active - BL Lacertae objects: general 



1 INTRODUCTION 

Blazars are a subclass of Active Galactic Nuclei (AGN) with 
a non-thermal emission covering most of the electromagnetic 
spectrum, i.e. from radio up to high-energy gamma-rays. 
Their broadband emission, that originates from a relativis- 
tic jet oriented close to the line of sight, is Doppler boosted 
and shows no evidence of spectral lines. The spectral en- 
ergy distribution (SED) of TeV-emitti ng blazars consists 
of two smooth, br oad components (e.g. lUlrich et al.l|l997l : 
iFossati et al.lll998l ). The first one extends from the radio up 
to the X-rays with a peak in the soft or hard X-rays, while 
the second one extends up to TeV energies, with a peak 
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energy around 0. 1 TeV, although this is not always clear 
(|Abdo et al.ll201ll ). 

Although the lower energy bump is attributed to the 
synchrotron radiation of relativistic electrons, the origin of 
the high-energy component is still under debate. Theoreti- 
cal models are divided into two classes: (i) leptonic and (ii) 
(lepto)hadronic, according to the type of particle responsible 
for the 7-ray emission. In the leptonic scenario, the high- 
energy component is the result of Compton scattering of 
electrons on a photon field; the seed photons can come either 
from synchrotron e mission of the same el e ctron population 
(SSC models) (e.g. iMaraschi erall Il992l : iKonopelko et all 
l2003n or fro m an 'externa l ' reg i on, i.e. from the ac- 
cretio n disk (JDermer et all Il992l : iDermer fc Schlickeiserl 
19931) or from the broad line region llSik ora et al.' 11994 : 



GhiseUini fc Madaul Il996l : iBoettcher fc Dermer 1998i ) (EC 



models). In the hadronic scenario, on the other hand. 
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the 7-ray emission is the resuh of proton or secondary 
lepton (th rough photohadronic interact i ons) synchrotron 
radiation IIMannheim fc BiermanrJ 1 19921 : lAhar onian 2000i; 
iMiicke et aLl2003h . For a recent review see lBoettchen (i2012f ) . 

The above models have been successfully applied for fit- 
ting the overall SED of TeV blazars by assuming stationary 
conditions, at least in most cases - this is particularly true 
for the class of hadronic models. However, one major feature 
of blazar emission is t he variability observ ed in almost all en- 
ergies (see for example lRaiteri et al.l2012l '). which further im- 
plies that stationary conditions might not apply. Moreover, 
variability can provide additional constraints on source mod- 
elling and therefore it can be used to lift the apparent degen- 
eracy between leptonic and hadronic models. As the above 
models use very different processes of widely varying cooling 
timescales, one expects that the system will react in a differ- 
ent way to variations on one or more source parameters, such 
as the injection rate of fresh particles. For example, one of 
the successes of SSC leptonic modelling is that it can repro- 
duce t he quadratic behaviour between the X-ray and TeV 
fluxes (|Mastichiadis fc Kirk|[l997l : iKrawczvnski et al.ll2002l ). 
The lack, thus far, of time-dependent hadronic models, has 
forbidden an analogous study of their expected radiative sig- 
natures. 

Recentlv lDimitrakoudis et al.l l|2012f ) have presented an 
one-zone hadronic model which is time-dependent. In the 
present paper we make use of the numerical code presented 
there in search for variability signatures in the context of 
these models. As an illustrative example, we have chosen 
to focus o n the relativel y recen t, high quality observations 
of Mrk421 iFossati et al] l|2008l '). that provide us with both 
spectral and temporal information. While we do not attempt 
to make a detailed spectral or temporal fit to the observa- 
tions, we use these as a springboard to examine the trends 
expected within the hadronic models. 

The present paper is organized as follows: In §2 we 
present the basic principles of the one-zone hadronic model 
and we comment on the possible options for fitting the SED 
of a blazar; the type of variations adopted in our simula- 
tions is also presented. In §3 we present our results starting 
with the multiwavelength (MW) fits obtained using three 
different (lepto) hadronic models. In subsections 3.2 - 3.3 we 
present the variability signatures obtained for each model 
caused by variations that we have applied on the injection 
compactness or on the maximum energy of particles respec- 
tively. We conclude in §4 with a summary and a discussion 
of our results. For the required transformations between the 
reference systems of the blazar and the observer, we have 
adopted a cosmology with iln-i — 0.3, Q,a = 0.7 and Ho = 70 
km s"^ Mpc"\ where the redshift of Mrk 421 z = 0.031 
corresponds to a luminosity distance Dl ~ 0.135 Gpc. 



2 THE MODEL 

2.1 General Principles 

In what f ollows we use the one - zone hadronic model as de- 
scribed in lDimitrakoudis et all l|2012l ') - henceforth DMPR. 
For completeness reasons we repeat here its basic points. We 
consider a spherical blob of radius R moving with a Doppler 
factor 5 with respect to us and containing a magnetic field 



of strength B. We further assume that ultra-relativistic pro- 
tons with a power law distribution of index pp between some 
energy limits 7p,min and 7p,max are injected into the source 
(we will be using Lorentz factors to denote proton or elec- 
tron energies throughout this paper). This injection can be 
characterised by a compactness 

i„j _ LpO-T 



'^ AnRnipC^ 



(1) 



where Lp is the proton injected luminosity and ctt is the 
Thomson cross section. 

Protons can lose energy via three channels: (a) syn- 
chrotron radiation, (b) photopair (Bethe-Heitler) and (c) 
photopion production. The effect that any of those three 
processes has on the proton distribution function depends 
on the specific parameters of the system, therefore all three 
have to be taken into account in a kinetic equation that 
also includes a proton injection and a proton escape term. 
Furthermore, since the above processes will create photons 
and other secondary particles (which will eventually decay 
to electrons and positrons, both of which we will hereafter 
refer to as electrons), one has to also follow the evolution 
of photons and electrons, by writing two additional kinetic 
equations for them q Assuming that the particles have a 
uniform distribution inside the source, one can write a sys- 
tem of partial, with respect to time and energy, integrodif- 
ferential equations, whose solution gives the corresponding 
particle distribution in addition to the multiwavelenrth pho- 
ton spectrum emerging from the source. 

The total number of free parameters used in this case is 
eight: The radius R of the source, the magnetic field strength 
B, the proton injection compactness ^p"^, the lower and up- 
per Lorentz factors of the proton power law distribution 
7p,min and 7p,max, as well as the proton index pp and the 
escape time from the source ip,osc. To these one should add 
the Doppler factor 5 of the emitting blob, which is used to 
boost the radiation and convert the source-frame quantities 
into observed flux. 

In the case where electrons are injected, in addition to 
protons, one could introduce the corresponding electron pa- 
rameters which are their injection compactness t^^ , which is 
related to their luminosity Le in the same way as described 
in eq. ((l| with rric replacing nip, the upper and lower cut- 
off of their injected spectrum 7c,min and 7o,max respectively, 
their slope pc and their escape timescale ic,csc. In order to 
reduce the number of free parameters one can assume that 

7p!i^iii — 7c,min — J- aUQ tp.csc — t^e,esc- 

Depending on the assumptions made about the time de- 
pendence of the parameters, the above scheme can be used 
to derive both steady-state and time-dependent solutions. 
Thus, if all parameters are constant in time the system will 
eventually reach a steady-state - note however that hadronic 
plasmas can become supercriti cal and in such cases they can 
exhibit limit cycle behaviour (jPetropoulou fc MastichiadisI 
l2012bl ) - hereafter PM12b. If, on the other hand, the sys- 
tem is in the subcritical regime and we allow for one or more 
parameters to have some explicit time dependence, then the 



^ Neutrons and neutrinos are also byproducts of photopionic in- 
teractions and one therefore should, in principle, write two more 
equations for them - see DMPR. For our present case , however, 
one can safely ignore them. 
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system will not reach a steady state but it will show continu- 
ous temporal variations, which will reflect the corresponding 
ones imposed on the input parameters and will have an im- 
pact on the produced spectrum. 

Therefore, one can first use the numerical code to obtain 
the SED of a source in a stationary state and then introduce 
perturbations in one or more of the fitting parameters to 
check the variability patterns in the MW spectrum. Similar 
methods have been applied in th e case of leptonic models b y 
iMastichiadis fc KirkI ^9^) and iKrawczvnski et all l|2002l ). 
However, in the framework of hadronic modelling we find 
that there are different combinations of radiative processes 
which can, in principle, give acceptable fits to the SED of 
blazars. Thus, it is possible that each fitting combination 
will have distinct temporal signatures which are worth in- 
vestigating. In what follows we make a qualitative discussion 
on the various options one has of fitting the SED of blazars 
using a hadronic model and then present our method re- 
garding temporal variations. 



2.2 Fitting the SED 

In contrast to the leptonic model where the relevant radia- 
tive processes are few, the hadronic model involves many 
processes which can make the radiated spectrum quite com- 
plicated. However, as it was shown in DMPR, there are cer- 
tain limiting, yet intuitive, cases, where the derived spec- 
trum has a particularly simple form. Thus, in the case where 
all features in the MW spectrum can be attributed to pro- 
tonqj, a case we will refer to as 'pure hadronic' or sim- 
ply 'model H', the spectrum has four distinctive features. 
One due to proton synchrotron radiation, two due to syn- 
chrotron radiation of electrons produced respectively in pho- 
topair and charged pion decay and one due to 7— rays from 
neutral pion decay. What is also interesting is that the fre- 
quencies where the three first peaks occur have a fixed ratio 
between them and scale as {m^/mp) : 1 : rj^^ where the 
value 77pe — 150 is deduced emp irically from the results of 
the SOPHIA code (|Miicke et al. 2000). So the frequencies 
of the first and third peaks are about eight orders of mag- 
nitude apart and therefore if the first one happens to be in 
the X-ray regime then, interestingly enough, the third one 
will be at TeV 7— rayo Clearly this coincidence needs some 
investigation as far as MW modelling is concerned. 

The introduction of a high luminosity leptonic compo- 
nent can change the picture significantly. If electrons have 
suitable parameters as to explain the X-rays, then there are 
two different options for fitting the TeV 7-rays with protons. 
The first is from the process described above, i.e. 7— rays are 
produced from a combination of the synchrotron radiation of 
electrons produced in charged pion decay and of the electro- 
magnetic cascade induced from neutral pion decay. We shall 
refer to this mechanism for TeV 7— ray production as 'pion 
induced' and we will denote the model with the acronym 



^ A primary leptonic component may also exist, provided its con- 
tribution to tlie synciirotron emission is much smaller than that 
of protons - see Appendix A for more details. 
^ Note that, as pointed by PM12b, the electromagnetic cascade 
that ensues from the absorption of the tt*^— decay 7— rays also 
contribute to the same energy regime. 



'LHtt'. As it turns out this model requires proton Lorentz 
factors not higher than 10^ and intermediate magnetic field 
strengths (~ 10 G). The same holds, in general, for model 
H. 

The second option is to fit the SED with very high pro- 
ton energies (7p,max ~ lO^'') and high magnetic fields (> 20 
G). In the latter case, proton synchrotron emission is re- 
sponsible for the TeV emission - this is the most popular 
fitting method of the relevant observations and we will refer 
to it as the 'proton-synchrotron' case or as 'LHs model'. Fi- 
nally, both LHtt and LHs models fall into the more general 
category of 'leptohadronic' cases. 

The above picture holds for low enough compactnesses. 
For higher ones photon-photon absorption, especially of the 
TT^-component that always peaks at very high frequencies, 
begins to dominate. However, as it was discussed in PM12b, 
the absorbed energy will be reemitted mostly at frequencies 
close to the peak produced by the radiation of charged pion- 
produced electrons, and only if the compactness becomes rel- 
atively high will the photon-photon absorption process be- 
come severe enough to distort the spectrum through electro- 
magmetic cascades and redistribution to lower frequencies. 
At still higher compactnesses the system undergoes a phase- 
transition and it becomes supercritic al: variou s radiative 
loops such as photopair-synchrotron ( Kirk fc Mastichiadia 
I992I) and photon quenching l|Stawarz fc KirkI 120071 : 



Petropoulou fc Mastichiadial201ll ) rapidly take energy away 



from the protons and redistribute it to electrons and ra- 
diation. In this case the system becomes very efficient, 
since a large fraction of the input proton luminosity is 
turned into radiation; however this is done at the cost 
of the spectral features described earlier which are de- 
stroyed by a very strong, non-linear photon-photon ab- 
sorption. Clearly one should seek successful spectral fits of 
MW blazar observations while the system is in the subcrit- 
ical regime; note howe ver that one can use supercr i ticali- 
ties in various contexts fMastichia dis fc Kazanaal2006l . l2009l : 
.Petropoulou fc Mastichiadia. 2012al ). 



2.3 An algorithm for inducing time variability 

DMPR have examined some examples of variability in the 
case of a pure hadronic model and showed that the system 
behaves like a SSC leptonic one, in the sense that proton 
synchrotron radiation produces the soft photons which act 
as targets for the photopair and photopion processes. Thus, 
in complete analogy to the leptonic SSC, hadrons interact 
with their own radiation. Consequently, one expects that 
variations in the proton injection parameter will produce a 
quadratic relation between the proton synchrotron and the 
pion induced components. 

While the inclusion of primary electrons introduces 
more free parameters in the model and facilitates the spec- 
tral fitting, it complicates the problem of inducing variability 
since now one should introduce two more free parameters, 
one of which relates the amplitude of variations between 
electrons and protons while the other relates their phases. 
In other words, while the pure hadronic model is expected 
to have a well-predicted temporal behaviour, since variations 
can be induced by changing a single parameter, the lepto- 
hadronic models will have, by neccessity, more complicated 
temporal patterns. 
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In order to simulate temporal variations we will adopt 
the following algorithm: 

(i) First we will obtain fits to a stationary/low state SED 
of a TeV blazar with each of the three aforementioned mod- 
els (H, LHtt, LHs). 

(ii) We will then introduce variations on some key param- 
eter (injection luminosity or maximum energy of particles) 
and follow in time the changes that these perturbations in- 
troduce to the SED of the source. 

(iii) Since most observations focus on the correlation be- 
tween X-rays and TeV 7-rays, we will focus on these energy 
bands and find the correlations expected in each of these 
three different scenarios. 

For concreteness, w e will apply this app roach to the ob- 
servations of Mrk 421 bv lFossati et al.l (120081 ). These sessions 
have produced a wealth of good quality, time-dependent 
data both at the X-ray (RXTE) and TeV (H.E.S.S.) regimes. 
However, we should warn the reader that we do not use 
the above data in order to make detailed fits but rather 
as a basis for comparing the general trends of our simula- 
tions with respect to them. Specifically, in the first step of 
the algorithm described above, we do not attempt a very 
detailed fitting of the SED but w e use the approach of 
IPetropoulou fc Mastichiadid (|2012al ) , where a family of pos- 
sible fits with low reduced x^ (below 1.5) can be considered 
as acceptable|j At any rate, it turns out that the expected 
temporal variations presented in §3 are largely independent 
of the value of the x^ of the particular fit. 

Motivated by the results of long-term variability 
studies of Mrk 421 and oth er prototype blazars (e.g. 
lEmmanoulopoulos et al.l (|2010l )). we have introduced, for 
the temporal variations, a random-walk type of change in 
one of the fitting parameters. Thus, we use 



tti+i = ai + (-1) 



int(5) 



i = 0,l,2. 



(2) 



where ^ is a uniformly distributed random number in the 
range (0,10). The integer parameter at is then scaled in such 
a way as to produce the following change in the parameter 
y we wish to vary 



yi = yo{l + fai), i = l,2, 



(3) 



where yi is the value of the parameter at time ti and / 
a multiplication factor - in all the examples shown in the 
present paper we have chosen / = 0.05, i.e. the parameter in 
question varies only by 5% between crossing times. For the 
values of y between two successive crossing times, we have 
chosen a linear interpolation scheme. Obviously no = 0, in 
order for the initial value of the parameter in question to 
match its corresponding value in the steady state fit yo- 



3 RESULTS 

We proceed next to show some characteristic results. We will 
begin by showing the spectral fits obtained and then we will 
focus on the variability induced on the models by varying 



* Clearly introducing more free parameters would have improved 
the fit but this would have acted against the scope of the present 
paper. 




E -10.5 




26 26.5 27 27.5 

logu (Hz) 

Figure 1. Top panel: Multiwavelength fit to the March 
22nd/23rd data of Mrk 421 using (a) a purely hadronic one-zone 
model (solid line), (b) a leptohadronic one-zone model where the 
TeV component of the SED is produced by synchrotron radiation 
of electrons resulting from photohadronic interactions (dashed 
line), and (c) as before but with the TeV component of the SED 
produced by proton synchrotron radiation (dotted line). For the 
parameter used see Table [T] . Bottom panel: The above zoomed 
in on the TeV energy range. 



first the injection compactness (i""^) and then the maximum 
energy (-y™^") of particles. 

3.1 Spectral fits 

Figure [l] (top panel) shows the spectral fits we obtained us- 
ing the three models described in the previous section to the 
pre-flaring state of Mrk 421 on the night of 22nd/23rd March 
2001. All data points used for the fitting are c ontemporane- 
ous; f or more details on the observations see iFossati et al.l 
(|2008r ). Different types of lines correspond to different mod- 
els as shown in the plot. A zoom in the TeV energy range 
is also shown in the bottom panel of the same figure. Note 
that, in all three models, the TeV 7-rays are fitted with the 
cutoff of the proton (LHs) or secondary lepton (H, LHtt) 
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Table 1. Parameters for the pure hadronic and leptohadronic fits 
to the pre-flaring state of Mrk 421 (see Fig. [T]|. 



Parameter symbol 


Model H 


Model LHtt 


Model LHs 


R (cm) 




3.2 X 10^^ 


3.2 X 10*5 


3.2 X 10^5 


B(G) 




20 


5 


50 


iJB (erg cm" 


-■') 


15.9 


1.0 


99.5 


<5 




16 


31 


21 


tt'r (hr) 




1.8 


0.9 


1.4 


'yp,niax 




8x 105 


4x 10^ 


4x 10^ 


Pp 




1.3 


1.5 


1.5 


*P 




1.6 X W-^ 


7.9 X 10"* 


1.6 X lO"'^ 


■7c,max 




- 


3 X 10* 


8 X 10^ 


Pe 




- 


0.7 


0.5 


C' 




- 


2 X 10"5 


5 X 10-5 


Up (erg cm" 


-3)1 


3.2 X 10"' 


1.6 X 10^ 


2.9 X IQ-i 


Uc (erg cm~ 


') 


2.3 X 10"* 


2 X 10"3 


3.4 X IQ--' 


Uj (erg cm" 


-'') 


1.2 


0.1 


3.6 X 10-1 


^j°t' (erg/s; 


\ 


6.9 X 10« 


1.3 X 10*8 


2.4 X 10** 



^ The listed particle and photon energy densities correspond to the 
pre-flaring fit, i.e. to a steady state of the system. 

synchrotron spectrum; in the latter cases (H and LHtt), 
synchrotron radiation from Bethe-Heitler electrons produce 
a distinctive broad spectral feature which lies between the 
proton synchrotron and the pion induced components. The 
fitting parameters of each model are listed in Table [l] In ad- 
dition, the comoving energy densities of protons (up), elec- 
trons (uc) and photons (uj) at the steady state are listed. 
The observed luminosity of the jet has been calculated ac- 
cording to 



P° 



TvR 3 /3c(llB + Up + Uc + U-y). 



(4) 



Although all three models can give us acceptable fits to the 
SED, the pure hadronic is the most energy demanding, hav- 
ing the largest ratio of particle to magnetic energy density 
(up/uB « 2 X 10^) and the highest jet luminosity. Clearly, as 
far as energy requirements go, the proton synchrotron model 
is, by far, the most economic of the three. 



3.2 Varying g^' 

One straightforward way of producing flaring activity is to 
vary the injection luminosity of the high-energy radiating 
particles. Thus, we start with the pure hadronic case and 
choose the proton injected compactness tp' to be the vary- 
ing parameter. Top panel of Fig. [2] depicts the X-ray (mid- 
dle) and the TeV lightcurves (bottom) centered at 10^* and 
2.5 X 10^® Hz respectively. We note that, in what follows, 
any reference to the TeV and X-ray fluxes will mean the 
integrated fluxes in the ranges (1.6 x lO'^^ — 10^^) Hz and 
(7.2 X lO^'^ — 3.6 X 10^*) Hz respectively. For comparison 
reasons, the injected proton luminosity (top curve) is also 
plotted after being rescaled. 

In the bottom panel of the same figure, the power spec- 
tral densities (PSD) of the injection and TeV light curves 
are shown with black and grey lines respectively. The corre- 
sponding PSD for the X-ray light curve is not shown, since 



it is similar to that of the TeVs apart from a normaliza- 
tion factor. While the PSD of the injection time-series can 
be fitted by a power-law with slope ~ —1.9, in agreement 
with the characteristic —2 slope of brownian noise, the PSD 
of the TeV light curve is best fitted with a broken power- 
law having the same slope with the injection's PSD at low 
frequencies (large timescales) but steepens significantly at 
higher frequencies (/ > 0.1). We note that same behaviour 
has been also detected in the two leptohadronic models. A 
study of the fitting parameters reveals that for protons pro- 
ducing the TeV 7— rays the condition icooi ^ icsc ~ icr 
holds, thus the system responds fast to the imposed vari- 
ations due to the small value of the proton escape time. 
On the other hand, electrons radiating in the X-rays have 
a much faster cooling than escape timescale. However, their 
PSD produces also a break at about the same frequency 
as in the TeV 7— rays case. This behaviour seems puzzling 
at first; one would expect that if the cooling of particles is 
fast, then the photon lightcurve (in our case the X-ray one) 
would track in detail the source variability. Moreover, one 
would expect to find a break in the PSD spectrum only for 
the TeV 7-rays, which is not the case. In fact, it is not the 
cooling timescale alone that determines the variability pat- 
tern of photons, but the minimum timescale between tosc 
and icooi. As is the case for all one-zone models, the fastest 
possible variation of the system is controlled by the photon 
crossing timescale icr- When min(tesc, tcooi) ~ tcr, as in our 
simulations, the photon lightcurves follow, in general, the 
source variations. However, the small timescale variations 
are smoothed out. Even in the extreme case of ultra fast 
cooling (icooi << tcr) photons cannot attain the full ampli- 
tude of particle variations - see Appendix B, and this results 
in a breaking frequency (which corresponds to a few icr) in 
the PSD. Another parameter that affects the shape of the 
output PSD is the amplitude of the imposed variations - 
see parameter / in eq. Q. For some very small value of /, 
which in our simulations is / — 0.001, we find that the lep- 
tohadronic system can follow the variations at injection, i.e. 
the PSD of the light curves has the same shape, apart from 
a normalization factor, with the one of injection. However, 
for so small f values the variability is at a very low level, 
and therefore for all practical purposes the source can be 
considered as not varying. 

For the two leptohadronic cases (Models LHtt and LHs 
in Fig. [T| we have implemented the same variatiorjj for tp^ , 
which was further extended to t^^ , since two types of parti- 
cles are being injected in the source. These simulations gen- 
erated light curves similar to those shown in Fig. [5] More 
information, however, can be deduced from flux-flux dia- 
grams, as the one shown in Fig. |3] where the TeV flux is 
plotted against the X-ray flux for all the cases discussed so 
far. First let us focus on model H (black solid line). While 
the X-rays follow linearly the proton synchrotron radiation, 
and therefore the proton luminosity, the TeV emission shows 
an almost quadratic dependence on tp' for the reasons ex- 
plained earlier - see also DMPR. Although at low fluxes one 
finds a clear quadratic dependence, this becomes fiatter at 
higher fiuxes due to the increasing effect of 77 absorption 



"" Not only the type of variation but also the series of random 
numbers used in all cases is the same, unless stated otherwise. 
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200 300 400 

t/tcr 
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Figure 2. Top panel: X-ray (dashed line) and TeV (dotted line) 
lightcurvcs for the pure hadronic fit shown in Fig. [T] obtained by 
varying tp' (thick solid line); the latter is shifted downwards by 7 
units in the y-axis for clarity reasons. For the adopted parameters, 
one tcr corresponds to ~ f .8 hrs in the observer's frame. Bottom 
panel: Power spectral density (PSD) of the injection time series 
(black line) and of the derived TeV light curve (grey line) . 



which gradually depletes the TeV regime. Another feature 
that is worth mentioning, is the tight correlation between the 
fluxes, which is derived with no requirements of fine tuning, 
same as in the SSC scenario. 

In general, the same trend holds also for the LHtt model 
(dashed line), while the LHs case (grey solid line) leads to 
a (sub)linear correlatiorjj. In the leptohadronic models we 
find no break in the TeV/X-ray correlations, since the effects 
of 77 absorption are less significant. This comes from the 

^ A linear regression fit to the data for the LHs model gives a 
slope ~ 0.7; we characterize this as a sublinear correlation. For 
a particular model, the exact slope of the Prev/^x correlation 
depends also on the energy bands considered. For example, for 
the LHs model we derive an exactly linear correlation between 
the peak fluxes of the two distinct emission components. 
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Figure 3. Plot of the TeV vs. X-ray fluxes obtained after varying 
only tp' (Model H) and both tp' and ^j."-" (Models LHtt & LHs). 
The fluxes are normalized with respect to their values of the pre- 
flaring state fit. Segments with different slopes are also plotted 
for reference. 



fact that, for both models, we have obtained good fits using 
much smaller tp' and higher Doppler factors than in the 
pure hadronic one - see Table[T] Finally, in all three models, 
no lags between the X-ray and TeV photons were detected. 

Previous analyses of Mrk 421 have noted the presenc e 
of lags in its spectra l evolution (e. g. Takahashi et al.l (|2000l ): 
iFossati et aP l|2008l ): ISingh et all (|2012h ). Such lags cannot 
be reproduced by one-particle population models, such as 
the purely hadronic one. However, in leptohadronic models 
where two primary particle populations are being injected 
into the source, these lags can be simulated by introduc- 
ing a shift N in the temporal variation of £^^' compared to 
^p"^, i.e. (£o"^)i = {tp')i+-N, where the subscript denotes that 
quantities are calculated at time fi (in units of ta). 

Figure |3] depicts the TeV vs. X-ray fiux obtained in 
the LHtt case. In both panels, the black line corresponds to 
A*' = 0, i.e. tp' and ^e"^ are completely correlated. The grey 
line in the top panel is the result of a simulation, where l]^' 
was obtained by a relative shift of 80 ta with respect to tp' . 
In this way we imposed less correlated variations on particle 
injection, having a Pearson's correlation coefficient r = 0.63. 
We have used a variety of positive correlated tp' , t^' by in- 
troducing various shifts, e.g N — 312, 512, 1024. In all cases, 
which we do not present here, we derived large correlation 
coefficients (r > 0.5) which show that the TeV and X-ray 
fiuxes retain their strong correlation. Then, we considered 
a more extreme case, where the proton and electron injec- 
tion have strong anticorrelation. For this, we have used two 
random number series with correlation coefficient r = —0.7. 
The resulting TeV/X-ray correlation is shown with grey line 
in the bottom panel. The respective diagram for the LHs 
model and A'' = 80 are shown in Fig. [5] 

In general, the introduction of a shift loosens the FtbV — 
Fx_ correlation. The effect, however, is more evident in the 
proton synchrotron case where the TeV/X-ray fiux correla- 
tion is almost destroyed even for A^ = 80 (r = —0.15). On 
the other hand, in the LHtt model, even if the input series 
are anticorrelated, we find a larger absolute value for the 
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Figure 4. Top panel: Plot of the TeV vs. X-ray fluxes obtained 
within the hHir model, after varying tp^ and t^^ with N = 
(black line) and N = 80 (grey line). Bottom panel: Same as above 
except for the grey line, which is obtained using as input for 
tp^ and i]^' two different anticorrelated time-series. Inset legends 
show the Pearson's correlation coefflcients for each case. 
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Figure 5. Plot of the TeV vs. X-ray fluxes obtained within the 
LHs model, after varying tp' and ^J."'' with N = (black line) 
and N = SO (grey line) . Inset legend same as in Fig. |4] 



correlation coefficient - see inset legend of Fig. |4l In other 
words, if the 7-rays are modelled by the emission of secon- 
daries produced by photohadronic processes and therefore 
any variations of FxeV reflect only indirectly the changes in 
the proton injection rate, the correlation between the fluxes 
in the X-rays and TeV 7-rays is partially retained. On the 
other hand, in the LHs model the variability observed in 
both X-rays and 7-rays reflects directly the variability pat- 
tern of the particle injection rate. Thus, if there is a degree 
of decorrelation in the injection it will be seen also in the 
flux-flux diagram. For the LHvr case (top panel) the intro- 
duction of a shift decreases the maximum/minimum flux 
values in both energy bands, since the 7-ray luminosity de- 
pends also on the number density of soft target photons, e.g. 
synchrotron photons in the X-rays. In the LHs case on the 
other hand, where the 7-rays are the product of proton syn- 
chrotron radiation, the 7-ray flux does not depend on the 
X-ray photons, as long as the X-ray photon number den- 
sity is low enough as not to cause signiflcant 77 absorption. 
Thus, the range of flux variations remains approximately the 
same. 



3.3 Varying 7max 

Observations of Mrk 421 during different periods of flaring 
activity indicate spectral evolution, and in particular spec- 
tral hardening i n the X-rays and/or in th e TeV 7-rays (e.g. 
pTakahashi et allbOOOl : iFossati et al.ll2008l ). Variations of the 
injection compactness alone cannot reproduce these observa- 
tional findings. This is exemplified in Fig. [51 where the spec- 
tral indices in the X-ray (grey line) and TeV energy bands 
(black line), are plotted against the corresponding fiuxes, 
for one of our fitting models (LHvr). The X-ray spectral in- 
dex remains approximately constant during flux variations, 
while the 7-ray part of the spectrum becomes softer during 
flaring events. This is to be expected, since the TeV observa- 
tions were fltted by the cutoff of the synchrotron component, 
in all three of our models. We note also that if 77 absorp- 
tion is signiflcant in the TeV regime, one expects to flnd an 
almost constant spectral index, even though the injection 
compactness varies. These two features, i.e. absence of spec- 
tral evolution in the X-rays and/or spectral softening in the 
7-ray regime, are also obtained for the pure hadronic and 
the proton synchrotron case. 

The second varying parameter that we have studied is 
the maximum energy of protons and electrons. We have kept 
constant the injection compactnesses and we have once again 
adopted the same random number series as in §3.2 for both 
7e,rrrax aud 7p,max. lu Figs. [7] and [5] wc show indicatively 
the results for the leptohadronic cases. For the proton syn- 
chrotron case, the linear flux-flux correlation is retained, in 
contrast to the photopion case. For the latter, we flnd a cor- 
relation between the TeV and X-ray fluxes, which is steeper 
than the quadratic we have obtained by varying £'^' and tp'. 
Actually, a flt to our results shows that -FreV oc -Fx'^- In the 
LHtt model the exact slope of the TeV/X-ray flux correla- 
tion is sensitive to the power-law exponent of the electron 
distribution. If the power-law spectrum of electrons at in- 
jection is flat, as in our LHtt fit (see Table[T|, any variations 



We define the spectral index as PV 
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Figure 6. Plot of the spectral index /3 as a function of the flux 
in the X-rays (grey line) and in the TeV energy regime (black 
line) for the photopion model, in the case where the injection 
compactness is the varying parameter. The fluxes are normalized 
with respect to their values obtained for the prc-fiaring fit shown 
in Fig. [T| 



Table 2. Slope of the TeV/X-ray correlations obtained for vari- 
ations on two model parameters. 



Parameter 


Model H 


Model LHtt 


Model LHs 


^inj 


linear to 
quadratic 


quadratic 


~ linear 


'Tmax 


quadratic 


quadratic to 
cubic 


~ linear 



of 7e,max result In non-negligible variations of the injection 
rate, since l]^^ is kept constant. Thus, the number density of 
synchrotron photons emitted by electrons with 7 < 7c, max, 
which serve as targets for protons, vary significantly. In such 
case, we expect that the 7-ray emission will show also large 
variations, since it depends on both the soft photon field 
and the proton distribution. On the other hand, when fitting 
the SED with a steeper power-law electron distribution (e.g. 
Pe = Pp = 1-5), the injection rate of lower energy electrons 
remains approximately constant. In this case, we obtain a 
quadratic relation between TeV and X-ray fluxes. 

Following the same procedure as in §3.2, we have then 
studied the effects of partially correlated variations of 7p,max 
and 7e,max on the TeV/X-ray correlation. We find that the 
degree of TeV/X-ray flux correlation is more easily destroyed 
in the LHs model. As for the LHtt model is concerned, we 
find, for the same shift, a more tight correlation in the case 
of variable injection than of variable energy cutoff. 

The TeV/X-ray correlations obtained so far (for cor- 
related input time-series) are summarized in Table (2] The 
symbol '~' is used to denote possible deviations from the 
exact trend due to effects that do not depend on the fitting 
model itself, such as 77 absorption and the particular choice 
of the energy bands. 

Figure [8] shows that an increasing flux in both X-rays 
(bottom panel) and 7-rays (top panel) is accompanied by 



Figure 7. Plot of the TeV vs. X-ray fluxes obtained after varying 
both 7p,max and 7e,max in Models LHvr & LHs shown in black and 
grey color respectively. The fluxes are normalized with respect to 
their values of the pre- flaring state fit. Segments with different 
slopes are also plotted for reference. 



a hardening of the spectrum; this trend is in good agree- 
ment with many observations regarding spectral variability 
of Mrk 421. Figure [8] reveals also a 'mirror' symmetry be- 
tween the photopion and proton synchrotron cases, i.e. the 
range of flux and spectral variations is larger in the TeV en- 
ergy band for the LHtt model and in the X-rays for the LHs. 
More specifically, the fact that Px varies less in the photo- 
pion case than in the proton synchrotron one can be used as 
diagnostic tool between the two models. The range of spec- 
tral variations in the two models is better displayed in Fig.|9l 
where the spectral index in the TeV energy range is plotted 
against the one in the X-rays. Although /3toV varies approx- 
imately between the same values in both models, spectral 
variations in the X-rays are more prominent in the proton 
synchrotron model. In particular, we find /?tcV oc 3.3/?x and 
/?TcV oc l.l/?x for LHtt and LHs respectively. 

The X-ray spectrum in the photopion case is relatively 
hard and shows smaller spectral variations due to the Bethe- 
Heitler component. If we artificially switch-off the channel 
of photopair production, find a fit to the pre-fiaring state 
of Mrk 421 and then produce variations to 7e,max/7p,max 
as before, we find that the range of variations for both /3x 
and Fx increases; specifically, /?x and logFx varry between 
(—1.7, —1) and (—9.9, —8.5) respectively, while their correla- 
tion law remains unaffected. Therefore, if the emission from 
the Bethe-Heitler process is neglected, no significant discrep- 
ancy between LHtt and LHs models is predicted. 



4 SUMMARY/DISCUSSION 

In the present paper we have examined the spectral and vari- 
ability signatures of the so-called (lepto)hardronic models. 
These models are routinely used to model MW observations 
of high energy blazars and constitute a viable alternative 
to the leptonic ones. According to standard practice, the 
low frequency part of the spectrum (usually up to the X-ray 
regime) is fitted by electron synchrotron radiation, while the 
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Figure 8. Plot of the spectral index /3 as a function of the fiux 
F in the TeV energy band (top panel) and in the X-rays (bottom 
panel). Black and grey lines correspond to the LHtt and LHs case 
respectively. For clarity reasons we have shifted the grey curve in 
the top panel by 0.2 units upwards. 
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Figure 9. Plot of the spectral index in the TeV energy range 
(/3xov) vs. the corresponding one in the X-rays (/3x) for the LHw 
(black line) and LHs (grey line) models. 



higher part is fitted by hadronic emission - most commonly 
proton synchrotron radiation. 

For the purposes of the present treatment we have used 
the time-dependent numerical code presented in DMPR. 
This models the photopair and photopion processes in great 
detail by using resu lts from the Monte Carl o code s of 
IProtheroe fc JohnsonI (|T996 ) and iMiicke et all (|2000l ) re- 
spectively for the production rates of secondaries. We have 
chosen, as an illustrative example, to fit the well mon- 
itored TeV blazar Mrk 421 duri ng a flaring state using 
the results of the 2001 campaign (jFossati et al.r 2008). For 
this, we followed the usual a lgorithm (for a simi l ar ap - 
proach to the leptonic cas e, see lMastichiadis fc KirkI (|l997h : 
iKrawczvnski et all ([2003)) where we first fitted the SED for 
a pre-flaring state and then we simulated variability by vary- 
ing some key parameter of the fit - in our case it was either 
the particle injection luminosity or the upper cutoff of the 
particle distribution. 

As explained in DMPR, hadronic models have three 
important processes which can lead to high energy pho- 
ton emission. These are proton synchrotron radiation, pho- 
topair, and photopion production. In principle all three 
could be used, in various combinations, in attempting a fit 
to the SED of a 7— ray blazar. In practice we found that 
there are three combinations which give satisfactory results 
to the Mrk 421 case: 

(i) The 'pure hadronic' (H) model: In this case there is 
no need for a leptonic component. The X-rays are produced 
from proton synchrotron radiation while the TeV 7-rays 
were 'pion induced' - with this we mean that 7-rays are 
produced from a combination of the synchrotron radiation of 
electrons produced in charged pion decay and of the electro- 
magnetic cascade induced from neutral pion decay 7— rays. 

(ii) The 'leptohadronic pion' (LHtt) model: The X-rays 
are produced from the synchrotron radiation of a primary 
leptonic component while the 7-rays are once again pion- 
induced. 

(iii) The 'leptohadronic synchrotron' (LHs) model: Here 
the X-rays are produced as in the previous model while the 
7-rays are produced by proton synchrotron radiation. 

All three models need high magnetic fields - as compared 
to the pure leptonic fits, and high values of 7p,max with the 
LHs model requiring the most extreme values of both pa- 
rameters. However, this is also the model that by far is the 
more 'economic' for the required energy density in relativis- 
tic protons - see Table [T] 

As a next step, we introduced random- walk type pertur- 
bations either in the compactness of the injected particles or 
in their upper cutoff. We have assumed that each consecutive 
value of these parameters is either increased or decreased by 
5% over its previous value. If many of these small amplitude 
perturbations pile up they can produce a 'statistical' fiare 
and we have chosen such an example to test the variabil- 
ity signatures of the models listed above (see Fig. [2| - note 
however that nowhere does the varying parameter exceed its 
initial value by more than a factor of 2. Feeding the values 
of the varying parameter to the code we produce simulated 
variations in the X-ray and TeV 7-ray bands and searched 
for possible correlations between their respective fiuxes and 
spectral indices. 

In the case of a varying injection compactness, the H 
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model gives us a clear quadratic dependence between the 
X-ray and TeV 7-rays - this can be easily explained from 
the fact that the H model resembles in many aspects to 
the SSC leptonic one (see DMPR). An interesting feature 
is that at high fluxes 77 absorption tends to linearize the 
correlation - see Fig. O The other two leptohadronic models 
need at least two more free parameters to fully specify their 
variability. One relates the amplitudes of variation between 
protons and electrons and the other their phases. Here we 
have examined the most simple cases, i.e. we have assumed 
equal amplitude changes between the two species which can 
be either contemporaneous or with a time shift of several 
fcr- In the no-lag case, the LHtt model produces a quadratic 
dependence between the fluxes, while the LHs a less than 
linear one. Note that because of the values of the fitting 
parameters (see Table [l]) 77 absorption is more severe in 
the latter case than in the former, a fact that explains the 
curvature seen in the model LHs curve. In the case where 
we have allowed time lags, then both models loosen up their 
correlation. We found however that in the case of the LHtt 
model a general quadratic trend remains while for the LHs 
model all correlation is lost. Finally, all three models seem 
to get steeper TeV spectra as the flux increases - see Fig[6l 

In the case where the maximum energy of the parti- 
cle distribution was varying, we found that the results are 
similar to the particle injection case, with the notable differ- 
ence that all models produce a spectral hardening to both 
X-rays and TeV 7— rays. We have also found that time shifts 
tend to decorrelate more the lightcurves, however the effect 
is stronger in the LHs model than in the LHtt. 

Obtaining the PSD of the light curves we find that they 
resemble the one of injection, up to a break frequency that 
corresponds to a few icr. In particular, for both TeV and X- 
ray light curves we find less power at the high-frequency part 
of their PSD (small timescale variations) when compared to 
that of the source. In other words, for the parameters used 
in the present work (see Table [l]), the photon field cannot 
react faster that a few ta to the imposed variations - see 
section 3.2 and Appendix B for more details. Note also that 
the break frequency of the PSD corresponds to ~ 0.5 — 1 
days, which does not c o ntradict other publish ed results (e.g. 
iTaka hashi ct al.' ('2000''):'Kataok a et al.l (120011 ) ). 

Although we have not addressed in detail the acceler- 
ation mechanism, the adopted modelling for the variations 
corresponds to a physical picture, where particles first are 
being accelerated close to a shock front (this region is a 
'black box' in our model) and then, they are being injected 
in the emitting volume. Thus, any changes that occur in 
the region where acceleration takes place, are later seen as 
changes in the characteristics of the particle injection mech- 
anism. These on their turn lead to an outbursting behaviour 
of the source. Studying the characteristics of this fiaring ac- 
tivity was the actual aim of the present work. In this re- 
spect, the introduced time-lags play the role of some in- 
trinsic difference in the acceleration mechanism between the 
two species which, however, retains some coherence. Two- 
zone model s whic h ha ve been applied to leptonic m odels 
(|Kirk et all l| 19981 ) and lMoraitis &: MastichiadisI l|201ll )) are 
more elaborate as far as particle acceleration is concerned, 
however they can become very cumbersome when treating 
hadronic processes. 

Concluding we can say that it is very interesting that 



all three models give very good x^ fits to the Mrk 421 data - 
the two leptohadronic models have slightly better fits, how- 
ever one should have in mind that they use many more free 
parameters. We note also that it is the first time that the 
pion-induced 7-ray component is used in fitting TeV data. 
In this case we found that a broad feature in the low to 
medium 7— rays is produced from the synchrotron radiation 
of Bethe-Heitler pairs (see Fig. 1). This feature is very re- 
strictive as far as fitting is concerned, but its presence might 
provide a decisive observational clue for the viability of this 
brand of hadronic models as it is not present neither in the 
LHs nor in the leptonic SSC models. We also found that 
proton synchrotron models are favoured as far as energetics 
is concerned. However, we note that the H and LHtt models 
produce in a much more natural way correlated variability 
between X-rays and TeV 7— rays which is mostly quadratic 
- that is they require equal amplitude variations between 
the particles and can allow for shifts in their phases. Proton 
synchrotron models, on the other hand, are inherently lin- 
ear and therefore, in order to produce quadratic correlations 
they need the rather special conditions that the electron in- 
jection varies quadratically with respect to the one of the 
protons. Furthermore, they require a tight phase correlation 
between the injection of the two species as all coherence is 
lost even in short time shifts. 

The analysis presented here deals with variability sig- 
natures in the framework of leptohadronic models, such as 
TeV/X-ray correlation, spectral evolution in the X-rays and 
very high energy (VHE) 7-rays and PSD of the correspond- 
ing light curves. These could be used as a diagnostic tool 
for differentiating between the models and therefore probing 
the nature of high-energy radiating particles. Observational 
requirements for this are the high temporal and spectral res- 
olution in the VHE part of the 7-ray spectrum along with 
contemporaneous X-ray observations. Great importance for 
the first requirement will be the future Cherenkov Telescope 
Array (CTA), since it will have high sensitivity and will 
provide us with large count rate light curves on a routine 
daily basis rath er than on exceptional fiaring events only 
(|Sol et al.ll2"oil ). 
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where 
dE 4 fm^Y 2 

= ^CTTCUB 7 • (A2) 

Assuming that protons at the injection have a power-law 
distibution, Up = Kp^~^ between 7p,min and 7p,max, with 
7p,min << 7p,max and p < 3, the integral above results in 

4 



ip*"" ~ ^(TTCUsKp 



2 3-p 
7p,max 



(A3) 



The normalization constant Kp can be expressed in terms 
of the proton injection compactness - see eq. ITJ in §2. In 
the case were proton cooling is not significant, the proton 
injection luminosity is approximately given by 



rinj 



where 



d7 7np(7) 



-K„F„ 



Fp 



Kp 



/p,max 



2-p 
' p,niin 



/p,niax 



2-p 



ifp<2. 



Thus, Kp is written as 



AnR 



2 /?inj 
P 



(JT-Fp 



(A4) 



(A5) 



(A6) 



For the electron distribution, on the other hand, we assume 
that all the injected power is radiated, i.e. L^^^ « ic"''- Thus, 
the ratio of the observed X-ray luminosities of the two pop- 
ulations is approximately given by 

1 









:40 



3-P /7p 



/mj 2 - 
*-p 



■P 



R ub 



106 1015 102 



.(A7) 



Assuming that y^ = y ~ 0.1, the above equation gives a 

rough estimation for t^^', which ensures that the emission 
features of primary electrons in the X-rays will be 'hidden' 
from the proton synchrotron component: 

, 2 - P X ^r 7p,max R Mb 



(2.5 X 10""): 



(A8) 



3-p 10-1 10-2 106 1015 102 
For the exact values listed in Table 1 we find tJ^' ~ 4 x 10 



APPENDIX B: SIMPLIFIED MODEL FOR THE 
SYSTEM'S RESPONSE TO VARIATIONS OF 
THE SOURCE 



APPENDIX A: LATENT PRIMARY LEPTONIC 
COMPONENT IN THE 'H' MODEL 

One of the models studied in the present work is the pure 
hadronic ('H'), where the two main emission components 
of the MW spectrum are attributed to protons. This as- 
sumption does not exclude the presence of a primary lep- 
tonic component, as long as its contribution to the total 
synchrotron emission in the X-ray regime is negligible. 

Thus, let us derive an approximate upper limit for the 
injection compactness of primary electrons £]^' . The proton 
synchrotron radiated power can be calculated by 



7-syn 
^P 



d7 np(7) ^ 



(Al) 



One of the factors that determines the response of a parti- 
cle/radiation system to the variations of the source is the re- 
lation between the variability timescale of the source, which 
in our simulations was taken equal to Iter, and the minimum 
of the cooling and escape timescales of particles. 

Here we show through a smiplified model, that even in 
the non-physical case of sub-dynamical cooling of particles, 
i.e. icooi << icr, the photon distribution cannot track ex- 
actly the variations of the source. For simplicity, we study 
the evolution of electrons and photons only. Let us assume 
that electrons, with number density n, are being injected 
into a volume with rate Q and escape in a typical timescale 



to 



while their characteristic cooling timescale is io,cooi 



p,rad 



energy losses of electrons are treated as catastrophic. Pho- 
tons, on the other hand, are escaping within t^r and are 
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being injected with a rate oc n/ic,cooi- Thus, the system is 
described by 



dn 
d7 



+ n{a + x) = Q(r) 



dn-,, 



+ n 



Jxn, 



(Bl) 
(B2) 



dr 

where r = t/tcr, a = tcr/io.osc, X ~ icr/ic,cooi and J is nor- 
mahzation constant. We model the electron injection func- 
tion as a series of pulses with constant duration {St = licr) 
and variable amplitude: 



where H{x) is the step function and 
Qi = Qj-2 +r, Qo = 2 and j = 2, 4, 6, ... 



(B3) 



(B4) 



and r is a uniformly distributed random number in the 
range [-3,2] while n = Ti_i + St and tq = 0. For initial 
conditions n(0) — n.y{0) — 0, the solution for electrons and 
photons is given by: 



Electrons: 



n = /i-i(r) 



n{Ti-i) + 



i-/i-i(^) 



a + X /i-l(T) 



i=l,3, 
i = 2,4. 



n(Ti-i 
where 

and each branch is valid in the time interval [ri_i, ri]. 
Photons: 



(B5) 



(B6) 



n^(ri_i)Si_i(r) + Jx 



a + X 



(l-iJi-i(r))+ (B7) 



JxGi_i(T) n(ri-i)- 



where 

Gi-i(r) 



(i5i-i(r)) 



a + X 



E,-i{r) 



1 -a- X 



and 



The solution for i = 2, 4, 6 is given below 

n(ri_i) 



n^(ri_i)Si_i(r) + Jx 



(B8) 
(B9) 

(BIO) 



Figure IBll shows the time-evolution of photons for three 
cases: (i) sub-dynamical cooling and fast particle escape with 
X = 100, a — 1 (top panel), (ii) slow cooling and fast escape 
with X ~ 0.1,0: = 1 (middle panel) and (iii) slow cooling 
and escape with x = 0.1, a = 0.1 (bottom panel). The nor- 
malization constant in examples (i)-(iii) was chosen to be 
J = 1, 20, 5 respectively. Note that even in the extreme case 
of ultra fast cooling (top panel) the photon lightcurve does 
not have exactly the shape of the injection pulses. Thus, the 
power spectral density (PSD) of the injection and photon 
time series will differ below some frequency /br. This effect 
will be even more prominent in examples (ii) and (iii), since 




Figure Bl. Response of the photon distribution to the vari- 
ations of electron injection for: x = 100, o = 1 (top panel), 
X = 0.1,0 = 1 (middle panel) and x = 0.1, a = 0.1 (bottom 
panel). The injection function is depicted with dashed Hues in all 



the photon lightcurve becomes smoother and loses the de- 
tails in small timescales - see middle and bottom panels in 
Fig. IBll In general, the temporal behaviour of the system 
including all physical processes, is satisfactorily described by 
one of the categories of the simplified model: 

• Fast cooling and escape, i.e. tcooi ^ iesc ~ icr- The 
leptonic component of our models belongs to this category, 
that corresponds to case (i) - top panel in Fig. IBll 

• Slow cooling and fast escape, i.e. icooi >> iesc = icr, 
which characterize the proton distribution in all our three 
models. Only in the LHs model, protons at the high-energy 
tail of the distribution have cooling timescales comparable 
to fcr- This category corresponds to case (ii) - middle panel 
in Fig. [BT] 

• Slow cooling and escape, i.e. fcooi ~ icsc 3> tcr- For the 
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fitting parameters adopted in the present worlc, neitlier the 
proton nor the leptonic distibution falls into this category, 
which corresponds to case (iii) - bottom panel in Fig. IBll 
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